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1. Introduction 

Stochastic processes in which independent particles scatter randomly at hnite speed 
and may occasionally propagate over a long distance in single bouts are known as Levy 
walks. Use of these models has become ubiquitous in the study of complex diffusive 
processes [1-5]. They are particularly relevant to situations such that the probability 
of a long jump decays slowly with its length [6]. The scale-free superdiffusive motion of 
Levy walkers [7] has been identihed as an efficient foraging strategy [8], spurring a large 
interest in these models, including with regards to human mobility patterns [9,10]. 
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A central quantity in this formalism is the distribution of free path lengths, 
which gives the probability that a particle propagates over a distance x between two 
successive scattering events. Its asymptotic scaling determines how the moments of the 
displacement asymptotically scale with time. Assuming the probability of a free path 
of length X scales asymptotically as one finds, in the so-called velocity picture 

of Levy walks, the following scaling laws for the mean squared displacement after time 
t [11-15]: 



0 

< 

a < 1, 

fVlogt, 

a 

= 

1, 


1 

< 

a <2, 

flogf, 

a 

= 

2, 

t, 

a 

> 

2; 


see below for a precise definition of this quantity. A scaling parameter value 0 < a < 2 is 
such that the variance of the distribution of free paths diverges, in which case the process 
is often called scale-free [9]. Correspondingly, the asymptotic divergence of the mean 
squared displacement gives rise to anomalous transport in the form of superdiffusion. 

In a recent paper [16], we considered Levy walks on lattices and generalized 
the standard description of the velocity picture of Levy walks, according to which 
a new jump event takes place as soon as the previous one is completed, to include 
an exponentially-distributed waiting time separating successive jumps. As emphasised 
in reference [16], this additional phase induces a differentiation between the states of 
particles which are in the process of completing a jump and those that are waiting to 
start a new one. We call the former states propagating and the latter scattering. The 
distinction between these states leads to a new theoretical framework of Levy walks in 
terms of multistate processes [17], whereby the generalized master equation approach to 
continuous time random walks [18,19] translates into a set of delay differential equations 
for the corresponding distributions. 

It is the purpose of this paper to show that a complete characterization of 
the solutions of such multistate processes can be obtained, which yields exact time- 
dependent analytic expressions of their mean squared displacement. These expressions 
can, on the one hand, be compared with the asymptotic solutions already reported 
in [16], and, on the other, also prove useful when the asymptotic regime is not reached, 
which is often the case with studies dealing with observational data. 

Such a situation arises when the probability of a transition from scattering to 
propagating states is small. A particle will then spend most of its time undergoing 
transitions among scattering states, performing a seemingly standard continuous-time 
random walk, only seldom undergoing a transition to a propagating state, during which 
it moves ballistically over a distance distributed according to the scaling parameter 
a. If a is small enough [a < 2), these occurrences, although rare, have a dramatic 
effect on the asymptotic scaling properties of the mean squared displacement, such that 
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a crossover from normal to anomalous diffusive transport is observed. This crossover 
time may in some situations, however, be much larger than the times accessible to 
measurements. Here we provide analytic results which give a precise characterisation of 
the emergence of an anomalous contribution out of a normally scaling one, showing how 
the asymptotic scaling (1) comes about. This is of particular relevance to the regime 
a = 2, for which the logarithmic divergence in time of the mean squared displacement 
is indeed very slow. Our results are tested and conhrmed by numerical simulations of 
the processes under consideration. 

The paper is organized as follows. In section 2, we introduce the multi-state 
description of Levy walks with respect to scattering and propagating states, and dehne 
the transition probabilities between them in terms of two parameters, one relating to 
the probability of a transition from a scattering to a propagating state and the other 
characterizing the asymptotic scaling of the free path distribution. The fraction of 
particles in a scattering state evolves according to a delay differential equation which is 
derived and solved in section 3. These results are exploited in section 4 to obtain the 
mean squared displacement of the processes. A comparison with numerical simulations 
is provided in section 5. Conclusions are drawn in section 6. 

2. Levy walks as multi-state processes 

We consider a continuous-time random walk on a square lattice which generalizes the 
standard model of Montroll & Weiss [18] in that it assumes displacements to non¬ 
neighbouring sites occur in a time span given by the ratio of the distance traveled to 
the walker’s speed n, which itself remains fixed throughout. In this sense, the model 
is similar to the so-called velocity picture of Levy walks [20] , the difference being that 
successive propagations, irrespective of their lengths, are separated by random waiting 
times, which we assume have exponential distributions. The model itself is not new and, 
in some sense, is a simplihcation of other models of Levy walks interrupted by rests; 
see reference [21]. As we show below, however, the combination of a discretized spatial 
structure and exponentially-distributed waiting times yields a novel description of the 
process in terms of delay differential equations amenable to analytic treatment. 

A natural distinction arises between the states of walkers which are moving across 
the lattice structure and those at rest. We call propagating the state of a particle which 
is undergoing a displacement phase and scattering the state of a particle at rest, waiting 
to start a new displacement. In the framework of intermittent random walks [22], the 
former state is usually referred to as ballistic and the latter as diffusive or reactive, 
depending on context. A scattering state may therefore be thought of as one associated 
with a local diffusive process bound to the scale ^ of the distance between neighbouring 
lattice sites. 

The interplay between the two states is as follows. A particle in a propagating 
state switches to a scattering state upon completing a displacement. A particle in 
a scattering state, on the other hand, can make transitions to both scattering and 
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propagating states; as soon as its randomly-distributed waiting time has elapsed, it 
moves on to a neighbouring site and, in doing so, may switch to a propagating state 
and continue its motion to the next site, or start anew in a scattering state. 

Consider a d-dimensional cubic lattice of individual cells n G The state of 

a walker at position n = (rii,... and time t can take on a countable number of 

different values, {k,j), specihed by an integer k E N, and direction j G {1,. .., z}, where 
z = 2d denotes the coordination number of the lattice, that is the number of different 
lattice directions. States (0,j) are associated with a scattering state, irrespective of 
direction j, while states (fc, j) with k>l refer to propagating states, with k being the 
remaining number of lattice sites the walker will travel in direction j to complete its 
displacement. 

Time evolution proceeds in steps, characterized by a waiting-time density function 
and a transition probability. A particle at position n in a scattering state will wait for 
a random time t, exponentially distributed with mean| tr, before updating its state 
to {k,j) with probability pk/z, simultaneously changing its location to n -|- e^, where 
Gj denotes the lattice vector (of length ^) associated with direction j. In contrast, a 
particle in a propagating state {k,j), k > 1, will change its state to (fc — l,j) after 
a time tb = i/v, simultaneously moving from site n to n -|- e^. The waiting times 
associated with scattering states are drawn from a standard Poisson process, while the 
renewal process generated by the combination of scattering and propagating states has 
arbitrary holding times, whose distribution is determined by the transition probabilities 
Pk- 

The waiting-time density of the process is thus the function 


^Jkit) 


T-^^e A; = 0, 

^d(A-tb), A; 7^ 0, 


( 2 ) 


where 5 d(-) denotes the Dirac delta function. When a step takes place, the transition 
probability to go from state (A;, j) to state {k\j') is 


'^{k,j),{k',j') 


z VfcG A; = 0, 

dk—i,k'djj' A; 7^ 0, 


(3) 


where 6.^. is the Kronecker symbol. 

A particle which makes a transition from the scattering state to a propagating state 
{k,j) will therefore travel a distance {k-\-l)i away in direction j, until it eventually comes 
back to the scattering state and may change directions at the next transition. 

We now introduce a characterization of the transition probabilities pk in terms of 
two parameters. The hrst, which we refer to as the scattering parameter, is denoted by 
0 < e < 1, and gives the total probability of a transition from the scattering state. 


f The subscript R in tr stands for residence as in “residence time.” In contrast, B in tb stands for 
ballistic as in “ballistic propagation time.” 
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'^Pk = e, (4) 

fc=i 

The remaining transition probability, 

Po = l-e, (5) 

is that of a transition from a scattering state into another scattering state. The value 
e = 0 thus corresponds to the absence of propagating states: the process is then a simple 
continuous-time random walk with transitions to nearest neighbouring sites only and 
exhibits normal diffusion with coefficient i‘^/{zT^). The opposite extreme, e = 1, assigns 
zero probability to transitions from scattering to scattering states, which means that 
every transition involves a displacement over a distance of at least two sites. Scattering 
states remain populated, however, due to the decay of propagating states. 

The second parameter, a > 0, is the scaling parameter of the transition probabilities 
Pfc, which controls their asymptotic behaviour, 

Pk oc k~°'~^ {k > 1), (6) 

and determines the scaling law of the mean squared displacement (1). The specific form 
of Pk has no effect on this scaling law, but is does affect the time-dependent properties of 
the mean squared displacement. To be specific, we consider in this paper the following 
double-telescopic form for the transition probabilities pk'- 

Pfc = e(l-2i-")-i[fc^-“-2(A: + l)^-“ + (A: + 2)i-“] (A: > 1); (7) 

its structure is particularly helpful for some of the computations presented below and 
is motivated by our study of anomalous transport in the infinite-horizon periodic 
Lorentz gas [23]. In this respect, the model is slightly different from that presented 
in reference [16], where a simple telescopic structure of the transition probabilities was 
used. 

For future reference, we define 


^k = ^p3, ( 8 ) 

j=k 

to be the probability of a transition to a state larger than or equal to A:, such that, in 
particular, uq = 1 and z/i = e, and, for the choice of transition probabilities (7), 


n = e(l - 2'-“)-' -(k+ 1)'““]. 


(9) 
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k 


= e{l - (1 - + I)'-” - (/; + 2)'-“]}, 


k 


( 10 ) 


=€(1 -2"-“)-'ll - (t+ 1)^-” + («: + 2)^-“ - 2(1:+ 2)"-“]. 


3. Fraction of particles in the scattering state 

A quantity which plays a central role in the analysis of the process generated by the 
waiting-time density (2) and transition probabilities (3) is the average return time to 
the scattering state, 


OO 



( 11 ) 


k=0 


The average return time is hnite only when§ a > 1. The process is then said to be 
positive recurrent. In the remainder, we restrict our attention to this range of parameter 
values. The null-recurrent case, when 0 < a < 1, for which the average return time to 
the scattering state diverges, is considered in reference [16]. 

The occupation probability of particles at site n and time t, P{n,t), is a sum of 
the probabilities over the different states, Pkj{n,t): 



( 12 ) 


j=l k=0 


In reference [16], we obtained the following set of delay differential equations for the 
time-evolution of these occupation probabilities: 



Z OO 


OO 



(13) 


k=0 



Z OO 


OO 



k'=0 


§ Generally speaking, pf. must decay asymptotically faster than k ^ for the average return time to be 
finite. 
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- crfc+fe/j(n - k'sj, t - {k' + l)rB) , 


( 14 ) 


where the inclusion of terms akj accounts for the possibility of external sources, such 
that f) > 0 is the rate of injection of particles at site n and time t in the state 

To simplify matters, we assume that these source terms are independent of the 
lattice direction and write crfcj(n,f) = z~^ak{'n,t). Furthermore, we will usually let 
crfcj(n, f) = 0 for k > 1 , which amounts to assuming that particles are injected only in 
a scattering state. The injection of particles in a propagating state is easily treatable 
as well, and will be considered explicitly in order to initiate the process in a stationary 
state of equations (13)-(14). We will, however, limit such considerations to this specific 
choice and avoid the possibility that source terms interfere with the asymptotic scaling 
of the process generated by particles initially injected in a scattering state, as might 
arise from alternative choices of injection rates of propagating states. 

Of particular interest is the fraction of particles in the scattering state. 


Z 

So{t) = (15) 

whose time-evolution is described by the following delay differential equation: 


OO OO 

So{t) = ^ PkSo{t - kTB) - er^^Soit) + ^ crfc(f - /ctb). (16) 

k=l fc=0 

Here, ak{t — kr-o) denotes a source term for the rate of injection of particles in state k, 
irrespective of their positions, that is, ak(t — kr^) = “ ^^b)- The first two 

terms on the right-hand side of this equation have the typical gain and loss structure 
of jump processes. Indeed, the second term corresponds to particles lost by scattering 
states due to transitions to propagating states, which occur at rate c/tr. Each such 
transition is gained back after a delay given by the length of the ballistic segment in a 
propagating state. The sum of those terms yields the first term on the right-hand side 
of equation (16). 

Before turning to the integration of this differential equation, we note that the 
asymptotic value of S'o(f), f —)■ oo, has a simple expression. By ergodicity of the process, 
the equilibrium ratio of particles in the scattering state is given by the ratio between 
the average time spent in the scattering state, tr, and the average return time to it, tq, 
equation (11), 


lim Soit) = —. (17) 

t —^CXD ^0 

Assuming a positive recurrent process, such as when the transition probabilities pk are 
specihed by equation (7) with a > 1, the average return time Tq = tr -|- e(l — 2^“")“^rB 
is finite, so that the asymptotic fraction of particles in a scattering state is strictly 
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positive. For the particular choice of the parameters pk, equation (7), the return time 
is given by equation (11) and equation (17) becomes 


lim S'o(t) = 

t^OO 


Tr 


rR + e(l-2i-«)-VB- 
Note, however, that the convergence to this asymptotic value follows a power law whose 
exponent goes to zero as a —)■ 1; see hgure la. Furthermore, So{t) converges to 0 when 
the scaling parameter falls into the null recurrent regime, 0 < a < 1. 

Correspondingly, the fraction of particles in the propagating states. 


tends to 


Z 

Sk{t) = ^Pfcj(n,f), 


(19) 


lim Sk{t) = -——j— , 

Tr + e(l — 


( 20 ) 


which follows from the observation that a particle that makes a transition to state k 
spends an equal amount of time in all states j such that 1 < j < k. 


3.1. Time-dependent fraction of particles in the scattering state 

The integration of equation (16) relies on the specihcation of initial conditions. For 
instance, the choice || 


So{t) = 


0, f < 0, 

1, t = 0, 

corresponds, when all particles are injected at the lattice origin, to the injection rate 


( 21 ) 


and will be henceforth referred to as the all-scattering initial condition. 

The time-dependent fraction of particles in the scattering state, S'o(t), t > 0, may 
then be obtained by the method of steps [24] , which consists of integrating the differential 
equation (16) successively over the intervals kr-Q <t<{k + 1)^6, A; G N, matching the 
solutions at the upper and lower endpoints of successive intervals. For the all-scattering 
initial condition (21), the fraction of particles in the scattering state is, for times t > 0, 

[i/'TBj k 

So{t) = a^n\k)rf[^it - kr^r, (23) 

k=l n=l 

where It/r^l (resp. [A/tb], used below) denotes the largest (resp. smallest) integer 
smaller (resp. larger) than or equal to A/tb, and each coefficient a(^n\k) is the sum of all 

II The discussion below can be easily generalized, by linear superposition, to processes where particles 
in a scattering state are continuously injected into the process. 
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possible combinations of products pi^... Pi„, where the sequences are such that 

*1 H- hin = k, 


Hri\k) = ^ 

l<2l ,...,2n<Al j = l 

ii-\ - \-in=k 

The contributions to a(^n\k) consist of all distinct ways of travelling a distance k inn steps, 
divided by their number of permutations; the derivation of equation (23) is provided 
in Appendix A. A comparison between the time-dependent fraction of particles in the 
scattering state (23) and the asymptotic value (17) is illustrated in figure la, where 
the difference between the two expressions is plotted vs. time for different values of 
the scaling parameter a and a specific choice of the scattering parameter e, with the 
transition probabilities pk specified by equation (7). 



(a) e(l - 2^-")-^ = 1 (b) e(l - 2^-")“^ = 5 x 10“2 


Figure 1: Convergence in time of the fraction of particles in a scattering state, 5*0 (f), 
computed from equation (23), to the asymptotic value (18). The transition probabilities 
are taken according to equation (7), with tb = 1 and tr = 1 -|- e(l — 2^““)“h The value 
of the scattering parameter is chosen such that e/(l — 2^“") = 1 (left panel) and 5 x 10“^ 
(right panel). The left panel shows the difference between So{t) and its asymptotic value. 
In both panels, the scaling parameter takes the values a = ?>/2 (green curve), 2 (cyan 
curve), and 5/2 (red curve). The right panel corresponds to a perturbative regime, such 
that So{t) can be approximated by a first-order polynomial in e (darker curves) whose 
first-order coefficient varies with time; see equation (29). The agreement between the 
exact and the approximate solutions improves as the value of e decreases. 


3.2. Small-parameter expansion 

Suppose that the scattering parameter is small, e 1, so that the overwhelming 
majority of transitions occur between scattering states, and only rarely does a particle 
make an excursion into a propagating state, which may, however, last long, depending 
on the value of the scaling parameter a > 1. 







Levy walks on lattices as multi-state processes 


10 


The first-order expansion of eqnation (23) yields the following approximation of S'q: 


So{t) ~ 1 - e-h > pk -, 

= S„([(/rBjrB)-t:ThlhB 




(25) 


k=\t/TB^ 


That is to say, for e <C 1, Soit) is a seqnence of straight line segments joining the valnes 
the fnnction takes at integer mnltiples of the propagation time tb, 


So{kTB) ^ 


L 


k_J2—ik-j) 

3=1 ^ 


(26) 


In this regime, S'o(f) converges asymptotically to the constant valne (17), which, for the 
choice of parameters (7), is given by eqnation (18), i.e., with all other parameters being 
hxed. 


lim So{t) ~ l-e(l-2^-“)-^ —. 

t—^OO 


For these parameters, we make nse of identities (10) to evalnate eqnation (26), 

So{kTB) ~ 1 - e(l - 2^ 
or, for general time valnes. 


"TR 


(27) 


(28) 


2'-“)-d- 

t'r 


([f/rej +1)^ " - {[t/TB\ +2) 


So{t) ~ 1 - e(l 

1 “ {[t/ 'TbJ + l) + ( [t/ tbJ 2) 

-2{[t/TB\+2) 


1—a 


JB 
tr l 


2-0 


, l — OL 


(29) 


A comparison between this approximate valne and the exact one, eqnation (23), is 
displayed in fignre lb, where the valne of e was taken to be large enongh that the cnrves 
remain distingnishable. 


4. Mean squared displacement 

Assuming initial injection of particles at the origin, the mean squared displacement of 
particles as a function of time, given by the second moment of the displacement vector 
r = £n, (r^)i = = P' XlneZ'* f), where = r • r and = n ■ n, evolves 

according to the differential equation 
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OCJ OCJ cxj 

— = T~^ X](2A: + l)i^kSo{t - /ctb) + ^(2/c - 1) crfc+fc'(0, t - fere), (30) 

k=0 k=l k'=0 

where, in contrast to the expression derived in reference [16], we have here added possible 
contributions from source terms (Tfcj(n, f), k > 1, which we assume to be concentrated 
at site n = 0 only. To conform with stationarity of the process, propagating states 
should be uniformly injected in the time interval —tb < f < 0 so they decay uniformly 
in the time interval 0 < t < tb- We thus let {k > 1) 

crA:(n, t) = Tg Vfc^n,o©(-t)0(t + ^-b), (31) 

where {fik}k>i is a sequence of positive numbers such that YlkLk = 1 and 0(.) is the 
Heaviside step function, such that the product 0(—f)0(f + tb) guarantees the uniform 
distribution of propagating source terms in the desired time interval. Equation (30) 
transforms to 

d ^ oo 

5Z(2fc + lWSa{t - /ctb) + rB^(2[t/rBj + l)^p\t/r^\+k- (32) 

fe=0 k=l 

Integrating over time, we obtain the mean squared displacement: 

LWb J «t/TB-k °° 

(n^)t = — V] (2fc + l)z/fc / dx So{xtb) +dx (2 [xj + l)pLa,j+fc. (33) 

k=0 k=l “^0 

There are two contributions to this expression; the first arises from the fraction of 
particles in a scattering state, and the second from particles initially distributed among 
propagating states. This term transforms to 


oo nt/TB °° L*/'^bJ 

dx (2[xJ + l)/xpj+fc = (2j ^)k'j+k 

k=l k=0 j=l 

oo 

+ 2{t/TB - [t/TB\){[t/TB\ + 1/2) . (34) 

k=l 


4-1- Stationary regimes 

The stationary initial condition^ i.e. such that the initial fractions of particles in 
scattering and propagating states are stationary solutions of equations (13)-(14), are 
realized if, in equation (31), we identify pk with the stationary fraction of particles in 
state k, 

_ _«4Tb_ 

tr+ e(l - 


( 35 ) 



Levy walks on lattices as multi-state processes 
cf. equation (20), and let 


12 


<^o(n, t) — p^oS]:,(t)6np, 




_Tr_ 

Tr + e(l - 2i-”)-'rB' 


( 36 ) 


With this choice, the mean squared displacement (33) yields the exact expression 


= 


t 


(1 - 21 “)rR + ctb 


Lh-T3j 


1 - 2^-“ + - (t - kTB){2k + l)[fc^-“ -{k + 1)^-“] 


k=l 


+ 


(1 - 2i-")rR + ctb 


I -^ + 2(VrB- LV^bJ) 

k=l 


[t/TBj + 1/2 

([t/TBj + 1)"-1 


(37) 


Treating n = ^/tb 3> 1 as an integer and letting Hn ^ = Ylk=i denote the generalized 
harmonic numbers, the two terms with sums over k in the above expression evaluate 
respectively to 


n 

^n-\n - k){2k + l)[k^-^ - {k + 1)^-“] ~ 1 - , (38) 

k=l 

for the term arising from the fraction of particles in a scattering state, and to 

n 

J2(2k - l)k^~° = , (39) 

k=l 

for the term arising from particles initially distributed among the propagating states. 
Depending on the value of the scaling parameter a, we have the following three 
asymptotic regimes. 

4.I.I. Normal diffusion For parameter values a > 2, the asymptotic properties 
of the harmonic numbers. 


lim ^ 

n^oo 

lim = C(a - 1). 

n^oo 

are such that only the terms arising from the fraction of particles in a scattering 
asymptotically contribute to equation (37): 



lim = — 

t 


■)l—a 


+ e[l -\- 2({a — 1)] 


(1 - 2i-“)rR + ctb 
which, up to a factor z~^, is the diffusion coefficient of the process. 


( 41 ) 
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4-1-2. Weak superdiffusion For the marginal parameter value a = 2, the 
harmonic numbers in equation (38) evaluate to 


//<"' = n, 

— logn + 7 + 

where 7 ~ 0.577216 is Euler’s constant. 

We therefore have the asymptotic limit of equation (37), 


(42) 


lim —— — — r{n^)t = -, (43) 

t^oo t \og{t/ tb) ' tr + 2erB 

which is due to the fraction of particles in a scattering state alone. This result 
is, however, of limited use because the asymptotic regime only emerges provided 
elog(f/rB) S> 1 (where 1 corresponds to the order of the sub-leading term), which may 
not be attainable, especially when the scattering parameter is small. For this reason, it 
is preferable to retain also the next-order term in equation (37), 


^ ^^2erB ^ 4e[log(t/rB) + 7 - 1/2 + 0(t/rB) ^] } , (44) 

where a fraction 4e/(rR -|- 2erB) is contributed by particles initially distributed among 
propagating states. 


4.1.3. Superdiffusion In the range of parameters 1 < a < 2, we substitute the 
scaling properties of the harmonic numbers. 


lim 

n^oo 

lim 

n—)-oo 

and obtain from equation (37) the limit 


(3-a) \ 
( 2 -«)-\ 


(45) 


lim 

t—^OO 


1 

(^/ub)^"" 



a — 1 1 1 2erB 

(2 — a)(3 — a) 3 — a (1 — 2^““)rR -|- ctb ’ 
_1 _2e7B_ 

(2 - a)(3 - a) (1 - + ctb ' 


(46) 


As emphasized by the hrst line of this equation, the leading-order coefficient is the sum 
of two distinct non-trivial contributions from the two terms on the right-hand side of 
equation (37): one from the fraction of particles in a scattering state and the other 
from particles initially injected in propagating states. As a consequence, the transport 
coefficient corresponding to the all-scattering initial condition differs by a factor a — 1 
from that above, corresponding to the stationary initial condition. This observation is 
consistent with the results of reference [25], where non-recurrent regimes of the scaling 
parameter a were also investigated. 
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The following exact expression of the mean squared displacement (33) is obtained after 
substitution of the fraction of particles in a scattering state, equation (23), and applies 
to the all-scattering initial condition (22): 


ib^Bj 




,-e(t-fcTB)/TR 
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j r 


E EE'^”“(”1-) T.r‘{r 
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ic-e-TB/m 


E E ‘’"“("I".) E |«‘ (?^ 

m=l n=l i=0 ■ ^ 
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([f/rej -k-my- (f/re - fc - 


(47) 


where, for consistency, one must interpret terms such as {j — m)° as unity, even when 
j = m. 



Figure 2; Growth in time of the rescaled mean squared displacement, {n‘^)t/t, computed 
from equation (47). The transition probabilities were set according to equations (7), 
with tb = 1 and tr = 1 -|- e(l — 2^““)“^. The value of the scattering parameter was 
chosen such that e/(l — 2^“") = 1 (left panel) and 5 x 10“^ (right panel). In both hgures, 
the scaling parameter takes the values a = 3/2 (green curve), 2 (cyan curve), and 5/2 
(red curve). For comparison, in the left panel, the asymptotic scaling values predicted 
by equations (41), (43) and (46) are displayed in dashed lines with matching colors. 
The right panel shows a comparison between the exact solution and an approximated 
solution, exact to first order in e (darker curves); see equation (53). 
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Figure 2 displays the results of explicit evaluations of equation (47) for three 
different values of the scaling parameter, a > 1, with e = 1 — According to 

the asymptotic results (41), (43) and (46), the scaling parameter value a = 2 gives 
rise to a logarithmic divergence of {n‘^)t/t, which separates the superdiffusive regime, 
for 1 < a < 2, from that of normal diffusion, for a > 2. A comparison between the 
exact solution and the asymptotic ones is displayed in hgure 2a. On the time scale of 
the hgure, convergence to the asymptotic regimes is convincingly observed only in the 
superdiffusive case. Note that while the computation of equation (47) up to t = SOtb 
takes less than one half hour of CPU time on a reasonably fast computer, doubling the 
time range would increase the required CPU time to over ten days. 

Although equation (47) gives the exact time-dependence of the mean squared 
displacement of the process for the initial condition (21), it does not give much insight 
into its time development. In particular, how to extract the large-time behaviour of 
the mean squared displacement and connect this result to the asymptotic scalings (41), 
(43) and (46) is not transparent. As rehected by hgure 2a, a numerical evaluation of 
equation (47) is necessarily limited to moderately large times, due to the number of 
terms involved. 

A regime of specihc interest which allows us to infer the emergence of the asymptotic 
scalings (41), (43) and (46) from the solution (47) is, however, that of small values of the 
scattering parameter e, i.e., such that transitions from scattering to propagating states 
occur only rarely. This is discussed below. We will otherwise have to resort to numerical 
simulations of the underlying stochastic processes to observe this convergence. Those 
results are presented in section 5. 

4-3. Small-parameter expansion 

For small values of the scattering parameter e, when transitions from scattering 
to propagating states are much rarer than transitions between scattering states, an 
expansion of equation (47) in powers of this parameter provides an approximate 
expression of the mean squared displacement from which the different asymptotic 
regimes discussed in section 4.1 can be inferred. In the anomalous regime of the scaling 
parameter, a > 2, terms constant in time, that carry a normal contribution to the mean 
squared displacement, may bring about contributions which, even for large times, may 
be much larger than that of terms diverging in time; an anomalous contribution to the 
mean squared displacement may then be masked by a normal one. 

To proceed, we substitute equation (28) into equation (33) and let t = kr-Q. Having 
dropped all terms of order and higher, we obtain 

{'n?)kT^ = — f da:|l-e(l-2^-")-^ —[1-([xj+1)^-"] 

'Tr Jo 


oo k oo 
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(48) 


where, in the last line, we have omitted terms that decay to zero as k grows large. 
Plugging into this expression the asymptotic forms of the generalized harmonic numbers 
(40), (42) and (45), 


C(a-l), a >2, 

hogA^ + 7-2, a = 2, (49) 

1 ( 2 - 54 — + 1 <“< 2 . 

we obtain approximations of the mean squared displacement (47), which can be 
compared with the asymptotic expressions (41), (43) and (46). 

In the regime of normal diffusion, a > 2, equation (48) reduces to the first order 
expansion in e of equation (41). For 1 < a < 2, the mean squared displacement (48) 
displays normal diffusion at short times and anomalous diffusion at large times. Thus 
in the weak superdiffusive regime, a = 2, equation (48) reduces to 


1 

t 
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l + 2e 


21og—+ 27 - 3 - — 
Tb tr. 


(50) 


which differs from the stationary expression (44) by a factor 4e/rR, subleading with 
respect to the logarithmically diverging term, identical in both expression. This 
difference stems from the absence of particles populating propagating states in our 
choice of initial condition. In the superdiffusive regime, 1 < a < 2, the divergence of 
the harmonic numbers, equation (49), dominates the large time behaviour. 
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(51) 


Here again we note a difference between this expression and equation (46), obtained in 
the stationary regime. At variance with the case a = 2, however, the difference occurs 
in the coefficient of the leading term on the right-hand side of equation (51), whose 
numerator is 1 — a instead of 1 in the stationary regime. The transport coefficient on 
the right-hand side of equation (51) therefore changes, depending on the choice of initial 
conditions. 

At large times, after the mean squared displacement crosses over from normal 
diffusion to superdiffusion, the above expressions are equivalent to 0(e) expansions of 
the asymptotic values (43) and (46). The value of the crossover time, tc, that separates 
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the short time normal diffusion from the large time anomalous diffusion can be inferred 
from the above expressions: 


fc ~ Tb X 


(2-a)(3-«) 
2(a-l) e 
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.exp , 


1 < a < 2, 

a = 2, 


(52) 


which can be large, in particular when a = 2. 

For short times, equation (48) can be improved by removing the assumption that t 
is an integer multiple of tb- One then finds 
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(53) 


A comparison between this approximate solution and the exact one, equation (47), is 
shown in figure 2b for different values of the scaling parameter, a. 


5. Numerical computations 

Numerical simulations of the process with rates (2) and probabilities (3) are based on a 
classic kinetic Monte Carlo algorithm [26] , taking into account the possibility of ballistic 
motion of particles in a propagating phase. 

A collection of independent walkers are initialized at time f = 0 at the origin of 
the two-dimensional square lattice 7?^ either in the scattering state, ko = 0, for the 
all-scattering initial condition, or in state fco ^ 0 with relative weights specified by 
equations (35)-(36), for the stationary initial condition. For each walker, we generate 
a sequence {(^n, jn), ^n}neN of successive states (fc„,_)„), G N, = 1,..., z = 4, and 
corresponding times as follows. 

When in a state kn-i = 0, the next transition is determined by drawing the following 
three random numbers. This is referred to as a scattering step. 

(i) The first random number, rj G [0,1], is drawn from a uniform distribution, and 
yields the state k^, such that Y!Xo Pa < V < YltloPa- 
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(ii) The corresponding waiting time, i.e., the time tn the particle waits in the state 
kn-i = 0 before the transition to state /c„ takes place, is obtained by drawing a 
second random number, whose distribution is exponential with mean tr. 

(iii) A third random integer, with uniform distribution among the set of lattice directions 
{1,...,4}, identihes the direction jn of the corresponding displacement by one 
lattice unit. 

Depending on the new state k^ another scattering step is taken if /c„ = 0 or, if 
/c„ > 1, a sequence of propagating steps jm), takes place, such that: 

(i) the state decreases by one unit, km = km-i — 1, 

(ii) the direction remains unchanged jm = jm-i-, and 

(iii) the time step takes value tm = tb, corresponding to the propagation time over a 
single lattice cell. 

The above loop repeats itself until m = n + k^ that is until km = 0, at which point 
another scattering step is taken. 

Keeping track of the positions of all walkers as functions of time, which we measure 
at regular intervals on a logarithmic timescale, measurements of the mean squared 
displacement {n‘^)t rescaled by 1/f are performed by taking averages of this quantity over 
the set of all walkers, which typically consists of 10® trajectories. Such measurements 
are reported below for the different scaling regimes in the positive recurrent range of 
the scaling parameter values, a > 1. 

Throughout this section, we set the transition probabilities according to equation 
(7) and, for convenience, change the scattering parameter e to 


1 - 2i-“ + e ■ 

We further let 



7-6 = 1 , 

2 (55) 

“ 2 - 5 ’ 

such that the asymptotic fraction of particles in the scattering state (18) becomes 
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t^OO 


2 + S' 


( 56 ) 
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With the choice of parameters (54)-(55), The mean squared displacement, for t large 
and with the all-scattering initial condition, is expected to scale as 


{n^)t 2 

t ~ 2 + 6 


1 -|- 5 (■(« — 1 ), 

l-F(5(logf-F7-2), 


1 + 5 


a—1 


(2-o)(3-o) 


t^ " + C(tt 



a > 2 , 
a = 2, 

1 < a < 2 , 


(all-scattering i. c.) 


(57) 

where, in the two anomalous cases, we kept terms constant in time to reflect the 
possibility that, when 6 is small, the normal term may not be negligible with respect 
to the anomalous one over a large range of times. This is particularly relevant for the 
marginal case, a = 2, where logf and 7 — 2 remain of comparable sizes throughout 
the time range accessible to numerical computations. In comparison, for the stationary 
initial condition, the above expression remains unchanged in the normal diffusive case, 
but is modihed in the two anomalous cases. 


{n^)i 


2 + 6 


1 + 5 C(« — 1), 

X < 1 + 5(logf + 7 - 1), 


1 + 5 


(2-a)(3-«) 


+ C(a - 1) 


a> 2, 

a = 2 , 

1 < a < 2 . 


(stationary i. c. 


(58) 


The hrst order expansion in 5 of the two anomalous regimes in equation (57) (in 
all-scattering initial condition) yields results equivalent to equations (50) and (51) 
respectively. 

A computation of the time evolution of the mean squared displacement in regimes 
of normal {a > 2) and anomalous (1 < a < 2) diffusion is shown in hgure 3, providing 
a comparison between the two initial conditions analyzed in section 4, with the two 
corresponding sets of asymptotic regimes given by equations (57) and (58). The scaling 
parameter values are set to a = 5/2 (hgure 3a), 2 (hgure 3b), and 3/2 (hgure 3c). The 
scattering parameter value is set to e = 1 throughout, such that transitions between 
scattering states are forbidden. As expected, the numerically computed mean squared 
displacement obtained from the stationary initial condition (magenta curves in hgures 
3a-3c) follow precisely the analytic result (37) in all three regimes. 

We further note that, in the regime of normal dihusion [a = 5/2), both data sets 
in hgure 3a display consistent convergence to the same asymptotic regime, given by the 
hrst lines of equations (57) and (58), limt^oo{'n‘^)t/t — 2.597. A similar result is observed 
in hgure 3b, where the convergence of the two data sets to the common leading value 
of equations (57) and (58), hmt_,.oo(n^)t/tlogf = 4/5, is apparent. The ehect of the 
dihering subleading terms, the one in (57) negative, the other in (58) positive (the 
latter value is about one half the absolute value of the former for 5 = 4/3), is, however. 
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(a) a = 5/2, 5 ~ 1.215 (e = 1) 
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(c) a = 3/2, 5 ~ 1.547 (e = 1) 


Figure 3: Numerical measurement of the time-evolution of the normally or anomalously 
rescaled mean squared displacement in the three distinct regimes of the scaling 
parameter. Time evolutions obtained from the stationary initial condition (magenta 
curves) are compared with those generated by the all-scattering initial condition (cyan 
curves). The dashed black curves correspond to the exact result (37) and the dotted 
lines to the asymptotic regimes (57) and (58), identical for regimes of normal diffusion, 
but otherwise differing according to the type of initial conditions. Insets: time evolution 
of the fraction of particles in the scattering state compared with the stationary value 
(56) (black dotted line). 
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manifest. Finally, in the superdiffusive regime, a = 3/2, figure 3c exhibits the two 
asymptotic regimes of the anomalously rescaled mean squared displacement given by 
equations (57) and (58), limt_,.oo(^^^)t/t^~“ — 0.582 (all-scattering initial condition) or 
1.163 (stationary initial condition), whose values differ by a 1 : 2 ratio for this value of 
the scaling parameter. 

5.2. Perturbative regimes of the scattering parameter 

The time evolution of the mean squared displacement in the perturbative regime of the 
scattering parameter, 5 1, is analyzed in hgure 4 for particles initially distributed 

in the scattering state at the origin, where the rescaled mean squared displacement 
is compared with the 0(e) solution, equation (53). Excellent agreement between 
the analytic and numerical results is observed throughout the time range when the 
parameter is small enough [6 = 10“^ hgures 4b, 4d, 4f). Numerical computations are 
also consistent with the asymptotic regime (57) for all cases. 

The scattering parameter value S = 10“^, shown in hgures 4a {a = 5/2), 4c (a = 2), 
and 4e {a = 3/2), is on the one hand small enough that, for short times, the computed 
mean squared displacement is barely distinguishable from the perturbative expansion 
(53). The effect of next order corrections is, on the other hand, apparent for larger times 
{t > 10^), beyond which a convergence to the asymptotic result (57) is observed. 

We note that the crossover time (52) is tc — 1.3 x 10^ for S = 10“^ (hgure 4c) 
and much larger yet for 6 = 10“^. For a = 2, subleading terms in equation (57) 
therefore dominate the logarithmically divergent term throughout the time range of 
measurements. Similarly, in the superdihusive case, a = 3/2, the respective crossover 
times (52) are p — 2x 10^ (hgure 4e) and p — 2.2 x 10"^ (hgure 4f), such that subleading 
terms in equation (57) remain important throughout the time range of measurements 
in both cases. 

We conclude this section by pointing out that the time range of numerical 
measurements such as presented in hgure 4 is limited by the hnite precision of the 
random number generator used to draw the transition probabilities pk and sets an 
ehective maximal scale of free hights, fcmax- This observation is analogous to the ehect 
of machine-dependent limitations recently reported in reference [27] and is particularly 
relevant to the regime of weak superdihusion. The ehective maximal scale fcmax induces a 
saturation of the logarithmic growth of the second moment for times larger than /cmax'^e- 
The process would thus become dihusive for times sufficiently large. Although this ehect 
can be pushed upward to larger times by increasing the precision of the random number 
generator, it cannot be eliminated altogether. 

6. Conclusions 

The inclusion of exponentially-distributed waiting times separating the successive jump 
events of LWy walkers leads to a natural distinction between propagating and scattering 
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Figure 4: Numerical computations (cyan curves) of the normally rescaled mean squared 
displacement in the three different regimes of the scaling parameter a for the all¬ 
scattering initial condition. The scattering parameter takes values in or near the 
perturbative regime: (a), (c), (e) 6 = 10-°, and (b), (d), (f) 6 = 10-°. The data 
are compared with the asymptotic solution (57) (black dotted curves) and the 0(e) 
solution (53) (dashed curves). Insets: time evolution of the fraction of particles in the 
scattering state with the stationary value (18) (dotted lines) and the the 0(e) solution 
(29) (dashed curves). 
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states, whose respective concentrations evolve in time according to a set of generalized 
master equations. In particular, the fraction of walkers in a scattering state obeys a 
linear delay differential equation with a countable hierarchy of delays whose analytic 
solutions were obtained. 

As opposed to classic methods based on the Fourier-Laplace transform of an integral 
kernel and the use of Tauberian theorems [11,13], the approach presented in this paper 
is based solely on the solutions of such delay differential equations. Our method thus 
yields a simple expression of the mean squared displacement of Levy walkers in terms 
of the distribution of free paths and the time integral of the fraction of particles in the 
scattering state. 

Both exact and asymptotic expressions of the mean squared displacement of walkers 
were obtained in regimes ranging from normal to superdiffusive subballistic transport. 
In these regimes, the mechanism through which a walker can change directions between 
successive propagation phases plays an important role in determining the values of 
the transport coefficients, whether normal or anomalous. The transition through a 
scattering phase brings about a description in terms of two parameters, the hrst 
specifying the typical timescale of scattering events as opposed to the timescale of 
propagation across an elementary cell, the second weighting the probability of transitions 
between scattering and propagating states. 

Relying on the stationary fractions of particles in scattering and propagating states, 
a detailed derivation of the transport coefficients, similar to that reported in [16], was 
given, exhibiting the precise effect of the two scattering parameters on these coefficients, 
as well as the influence of the initial distribution of walkers [25]. Our formalism also 
yields the exact time evolution of the mean squared displacement, which was investigated 
when particles are initially found in a scattering state. 

The comparison between the exact and asymptotic solutions is generally not 
immediate, but is fairly straightforward in perturbative regimes such that the likelihood 
of a transition from scattering to propagating states is small. A case in point - though 
by far not the only application of the present results - is given by the inhnite-horizon 
Lorentz gas [28], for which the scaling parameter of the distribution of free paths takes 
on the marginal value a = 2. As one might expect, measuring the logarithmic divergence 
in time of the rescaled mean squared displacement is typically hindered by dominant 
constant terms. This is particularly so in the regime of narrow corridors, which yields 
a stochastic description in terms of multistate Levy walks such that the scattering 
parameter is small, e -C 1 [23]. The parameter e induces a further separation between 
the time scales of the scattering and propagating phases, tr S> tr. To accurately 
model the scattering phase is therefore of paramount importance to understanding the 
dynamics of the inhnite-horizon Lorentz gas in terms of a Levy walk. 

The results obtained in this paper provide the framework to transposing such results 
to a range of applications exhibiting other regimes of transport, such as have been 
studied in the context of optimal intermittent search strategies [22]. Our theoretical 
framework brings about new perspectives to extend the study of intermittent walks to 
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power-law distributed ballistic phases, which may be relevant to the increasing body of 
literature on optimal search strategies [8,29,30]. 
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Appendix A. Time-dependent fraction of particles in the scattering state 

To derive equation (23), we note that a(i|fc) = pk and, for 2 < n < k, 
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By recursive application of this formula, we obtain 
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which is equivalent to equation (24). 

Taking the derivative of equation (23), we verify equation (16): 
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We have to show that the last term on the RHS vanishes, which amounts to proving the 
identity 
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= 0, (A.4) 


or, after rearrangement. 
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The second term in this expression, by successively changing the index j to j — k, then 
swapping the sums over k and j and exchanging the indices j and k, transforms to 
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Now swapping the sums over j and n in the last line and plugging this expression back 
into equation (A.5), we obtain, after factorization of the common factors, the condition 
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(A.7) 


which holds by identity (A.l), thus completing the proof that S'o(f) as specified by 
equation (23) is the solution to equation (16) for the initial condition (21). 
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